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Abstract 

Specific heats of quantum systems with symmetric and asymmetric double- well potentials have 
been calculated. In numerical calculations of their specific heats, we have adopted the combined 
method which takes into account not only eigenvalues of e n for < n < N m obtained by the 
energy-matrix diagonalization but also their extrapolated ones for N m + 1 < n < oo (N m = 20 or 
30). Calculated specific heats are shown to be rather different from counterparts of a harmonic 
oscillator. In particular, specific heats of symmetric double-well systems at very low temperatures 
have the Schottky-type anomaly, which is rooted to a small energy gap in low-lying two-level 
eigenstates induced by a tunneling through the potential barrier. The Schottky-type anomaly is 
removed when an asymmetry is introduced into the double-well potential. It has been pointed out 
that the specific-heat calculation of a double-well system reported by Feranchuk, Ulyanenkov and 
Kuz'min [Chem. Phys. 157, 61 (1991)] is misleading because the zeroth-order operator method 
they adopted neglects crucially important off-diagonal contributions. 
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I. INTRODUCTION 



Double-well (DW) potential models have been employed in a wide range of fields in- 
cluding physics, chemistry and biology (for a recent review on DW systems, see Ref. PQ). 
We may classify quantum DW models into three categories: exactly solvable, quasi-exactly 
solvable, and approximately solvable models [2]. In the exactly solvable model, we can de- 
termine the whole spectrum analytically by a finite number of algebraic steps. In contrast, 
we can determine a part of the whole spectrum in the quasi-exactly solvable model. In other 
models, eigenvalues are obtainable only by approximate, analytical or numerical method. 
Examples of exactly solvable models include the double square-well potential and the Man- 
ning potential [3]. The Razavy potential [U E] expressed by hyperbolic functions belongs 
to the quasi-exactly solvable models. In this paper, we pay our attention to two types of 
approximately solvable models with a quartic potential [6H9] and a quadratic potential per- 
turbed by a Gaussian barrier [101 [TT] , which are hereafter referred to as model A and model 
B, respectively. These models have been commonly adopted for studies of tunneling and 
stochastic resonance in DW systems. It is, however, curious that studies on their thermo- 
dynamical properties are scanty [T2HT5] . Feymann and Kleinert applied the path-integral 
method to a calculation of an effective classical partition function of DW systems [121 H3] • 
Okopihska studied the effective potential, employing the optimized and mean-field expan- 
sions of a path-integral representation for the partition function [13]. The specific heat of 
a DW system was calculated by Feranchuk, Ulyanenkov and Kuz'min (FUK) JT5] with the 
use of the zeroth-order operator method (ZOM) [16], related discussion being given in Sec. 
IV. 

It is the purpose of the present paper to study the specific heat of quantum systems 
with symmetric and asymmetric DW potentials. Various kinds of analytical and numerical 
methods for approximately solvable DW models have been proposed to evaluate their eigen- 
values [6TfT2"l ITIHTB] . In the present study, we evaluate them by a numerical diagonalization 
of the energy matrix with a finite size N m (= 20 and 30). Eigenvalues e n for < n < N m 
are sufficient for a study of thermal properties of DW systems at very low temperature 
near T = K. However, they are insufficient for describing thermodynamical properties at 
elevated temperatures, as explicitly shown shortly (Figs. [3] and [9]). In order to overcome 
this deficit, we adopt the combined method in which we include additional eigenvalues e' n , 
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extrapolating to a larger n (N m + 1 < n < oo) such that they lead to results consistent with 
classical statistical calculations. Taking into account both extrapolated eigenvalues as well 
as those obtained by the energy-matrix diagonalization, we may obtain reasonable specific 
heats at both low and high temperatures. 

The paper is organized as follows. In Sec. II, we will describe the adopted combined 
method for the energy-matrix diagonalization and extrapolated eigenvalues. In Sec. Ill, the 
combined method is applied to DW systems with a quartic potential (model A), a quadratic 
potential with Gaussian barrier (model B) and the DW potential adopted by FUK (FUK 
model) [15J. In Sec. IV we critically examine a validity of the specific heat calculated by 
FUK [15] with the use of ZOM [16]. Specific heats of the triple- well system are studied also. 
Sec. V is devoted to our conclusion. 



II. THE COMBINED METHOD 



A. Classical statistical calculation 



We consider a system whose Hamiltonian is given by 

+ (1) 

where m and U(x) are a mass of a particle and a DW potential, respectively. The classical 
partition function is given by 

with 

ZM = / e-^W dx, (3) 



oo 



where /3 (= l/Zc^T) denotes an inverse temperature. From the calculated partition function 
Z({3), we obtain the specific heat C and entropy S by 

c - fr 

S=^(E-F), (5) 
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where 



F = ~lnZ03). (7) 



B. Quantum statistical calculation 

In order to make a quantum statistical calculation, it is necessary to evaluate eigenvalues 
of a given Hamiltonian of Eq. (JT]), for which the Shrodinger equation is given by 

h 2 d 2 



+ U{x) 



V(x) = EV(x), (8) 



2m dx 2 

$?(x) and E standing for eigenfunction and eigenvalue, respectively Various approximate 
analytical and numerical methods have been proposed to solve the Schrodinger equation 
[6-llJ. We evaluate eigenvalues, treating the Hamiltonian H as 

H = H + U(x) - U (x) =H Q + V{x), (9) 

with 

2 

#0 = £- + U (x), (10) 
2m 

V(x) = U(x) - U (x), (11) 
%(«) - (12) 

where uq denotes the frequency of a harmonic oscillator. Eigenfunction and eigenvalue for 
H are given by 

1 fmu y/4 ( mu x 2 \ ( Imuo \ 

E 0n = (n + hw , (14) 

where H n (x) stands for the Hermite polynomials. Expanding the eigenfunction ty(x) in 
terms of (j) n (x) 

oo 

= ^c n n (x), (15) 

?1=0 
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we obtain the secular equation for {c„} expressed by 

oo 

E °m = ^ H mn C n , (16) 



n=0 

with 



/oo 
(j) m (x)V(x)(j) n (x) dx, (17) 
■oo 

where c n denotes an expansion coefficient. The method mentioned above is not new, and 
equivalent or similar ones have been adopted in Refs. [6H9]. It is possible to calculate H mn 



by MATHEMATICA [T7]. We may diagonalize Eqs. (16) to obtain eigenvalues e n (n = to 



N m ), and evaluate the quantum partition function given by 

Z((3) = Tr e~ pH = ^V^ 6 ", (18) 

n=0 

where N m stands for the maximum eigenvalue. 

As will be shown shortly (Figs. [3] and [9]), eigenvalues for < n < N m (N m ~ 20 — 30) 
are sufficient for a study of low-temperature thermodynamical quantities, but insufficient 
for high-temperature ones. We adopt the combined method in which we include not only 
eigenvalues e n for < n < N m obtained by energy-matrix diagonalization but also their 
extrapolated ones for n > N m given by 

e' n = A n r hw for N m + 1 < n < oo, (19) 

with parameters A and r. An exponent r is chosen such that the high-temperature specific 
heat calculated with combined eigenvalues is consistent with the classical specific heat as 



follows. A simple calculation of the partition function with Eq. ( 19 ) in the high-temperature 
limit leads to 



Z(f3) J exp[-f3Az r ] dz = ((3A)- l,r V (l + 



(20) 



yielding the specific heat 



C = ( ' ) (21) 



which should be in agreement with the classical specific heat obtained by Eqs. Q and (|6]). A 



prefactor A in Eq. (19) is chosen such that e' for n > N m + 1 becomes a good extrapolation 



of eigenvalues e n for < n < N m evaluated by the energy-matrix diagonalization. Then the 
resultant quantum partition function is given by 

n=0 n =N m +l 

where the first and second terms express contributions from eigenvalues {e n } derived by the 
energy-matrix diagonalization and from extrapolated eigenvalues {e' n }, respectively. 

III. APPLICATIONS TO MODEL POTENTIALS 
A. A quartic DW potential (model A) 

1. The symmetric case 

First we apply our combined method to model A with the symmetric DW potential given 



U{x) = '^J^~ 4f> (23) 



by 

which has stable minima at x = ±xo and an unstable maximum at x = 0. The height of 
the potential barrier is A = U(0) — U(±x ) with U(0) = mujQxl/8 and U(±x ) = 0. For 
a later purpose of an energy- matrix calculation, we have chosen U"(±Xq) = mw^ such that 
the potential U(x) has the same curvatures at the minima as the harmonic potential Uq(x) 



given by Eq. (12). For numerical calculations, we assume m = 1.0, uq = 1.0 and xo = 3.5 
[T7j . for which If' (xq) = U' Q (0) = 1.0 and A = 1.531. The symmetric DW potential in 
model A is plotted by the solid curve in Fig. [TJ where the harmonic potential given by Eq. 



(12) is shown by the dashed curve: chain and double-chain curves will be explained later 
[Eq. j2§]. 

We have calculated the classical specific heat which is shown by the chain curve in Fig. 
[2| The calculated specific heat in the high-temperature limit of T — » oo reduces to 

C = (\ + \) k B = \ k B , (24) 

because we obtain C{T)/k B = 0.690, 0.7293, 0.7434 and 0.7473 for k B T/hw = 10.0, 100.0 



and 1000.0, respectively. The first (1/2) and second (1/4) terms in Eq. (24) express con 



tributions from momentum (p) and coordinate (x), respectively, the latter being due to the 



lOr 




FIG. 1: (Color online) The symmetric quartic potential [Eq. ( |23[ )] (solid curve) and asymmetric 
ones [Eq. (|28|)] with d = 0.01 (chain curve) and d = —0.01 (double chain curve) of model A, the 



dashed curve expressing the harmonic potential [Eq. (12)]. 



quartic power of the potential. Indeed, in a system with a quartic potential of U(x) = x 4 /4, 
the coordinate contribution to the classical specific heat becomes (1/4) ks (the Virial theo- 
rem). 

For a special case of the DW potential 



U(x) 



x 

T 



we obtain the analytical expression for Z x ((3), 



z x {P) 



71 



J/8 



I i 



X 

~2 



P 



(25) 



+ h 



P 



(26) 



which yields C{T)/k B = 0.7236, 0.7416, 0.7473 for T = 10, 100 and 1000, respectively, I n {z) 
denoting the modified Bessel function of the first kind. 
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FIG. 2: (Color online) Temperature dependences of classical (chain curve) and quantum specific 



heats (solid curve) of model A with the symmetric DW potential [Eq. ( 23 )] , dashed curve expressing 
quantum specific heat of a harmonic oscillator (HO). The inset shows an enlarged plot of the 
quantum specific heat at very low temperatures with the Schottky-type anomaly. 



Matrix elements H mn of Eq. (17) are finite for pairs of \m — n\ = 0, 2 and 4. Figure 
[3]^a) shows eigenvalues {e n } obtained for N m = 20 (open circles) and 30 (filled circles). 
Eigenvalues for n < 20 are almost the same for N m = 20 and 30. Calculated eigenvalues e n 
for n = 0, 1, 2 and 3 are 0.476188, 0.478131, 1.2695 and 1.3514, respectively, with N m = 30. 
Eigenvalues of eo = 0.476188 and e\ = 0.478131 originate from two eigenvalues of Eqq = 0.5 



in Eq. (14) at two minima with symmetric and antisymmetric wavef unctions. They are 
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FIG. 3: (Color online) (a) The n dependence of eigenvalues e n of model A with the symmetric 



DW potential [Eq. (23)] obtained by the energy-matrix diagonalization for < n < N m (circles) 
and by an extrapolation given by e' n = 0.25 n 4//3 Huiq for n > N m (dashed curve), the inset 
showing eigenvalues obtained by the energy-matrix diagonalization with N m = 20 (open circles) 
and N m = 30 (filled circles), (b) The temperature dependence of the quantum specific heat 
calculated with eigenvalues e n for < n < N m with N m = 20 (chain curve) and N m = 30 (dashed 
curve) , the solid curve expressing the result with combined eigenvalues shown by the dashed curve 
in (a). 

quasi-degenerated states with a small gap given by 5 = ei — e = 0.001943, which is induced 
by the tunneling effect through the potential barrier. Similarly, eigenvalues for n = 2 and 3 
are also quasi-degenerated as given by e 3 — e 2 = 0.0819. 

Specific heats calculated with the use of these eigenvalues e n for < n < N m with 
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FIG. 4: (Color online) Temperature dependences of classical (chain curve) and quantum entropies 



(solid curve) of model A with the symmetric potential [Eq. (23)], and the quantum entropy of 
a harmonic oscillator (HO) (dashed curve), the inset showing an enlarged plot of the quantum 
entropy at very low temperatures. 



N m = 20 and 30 are plotted by dashed and chain curves, respectively, in Fig. |3[b). They 
are in good agreement at /c#T /hu < 0.5 but significantly different at /c B T '/hu > 1.0. This 
implies that eigenvalues for < n < N m (= 20 and 30) are insufficient for a study of the 
specific heat at elevated temperatures of ksl '/huo > 1.0. 

We adopt the combined method with extrapolated eigenvalues given by 

e n = 0.25 n 4/3 hu for N m + 1 < n < oo, (27) 
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where an exponent of r = 4/3 chosen by C/ks = 1/r = 3/4 in Eqs. (21) and (24) is 
consistent with the WKB type analysis for a large n. Extrapolated eigenvalues given by 
Eq. (27) are plotted in Fig. [3]^a). The quantum specific heat of model A calculated with 



combined eigenvalues is shown by the solid curve in Fig. [2] [or Fig. [3[b)] [18J. The quantum 
specific heat is rather different from the classical one at low temperatures as expected. A 
closer inspection of the quantum specific heat reveals that C(T) has an anomalous peak at 
very low temperature at kgT '/ftiuo ~ 0.001 ~ 5/2, as shown in the inset of Fig. [2] It is the 
Schottky-type specific heat arising from low-lying two-level eigenvalues of eo and ex whose 
energy gap is induced by a mixing through a tunneling. Although quantum and classical 
specific heats do not well agree at ksT/heo ~ 10 in Fig. |2j both reduce to (3/4) &b in the 
high-temperature limit 

For a comparison, we show the quantum specific heat of a harmonic oscillator by the 
dashed curve in Fig. [2j The quantum specific heat of model A is not dissimilar to that of a 
harmonic oscillator at low temperatures except for the Schottky-type anomaly. However, the 
high-temperature specific heat of model A given by C = (3/4) fc^ is different from Cho — 
of a harmonic oscillator. 

Temperature dependences of classical and quantum entropies of model A are shown by 
chain and solid curves, respectively, in Fig. [4] where the quantum entropy of a harmonic 
oscillator is plotted by the dashed curve. With decreasing the temperature, the quantum 
entropy decreases but seems to remain at 0.69 (~ In 2). The inset of Fig. [1] shows that it 
furthermore decreases below k%T '/ hwo ~ 0.002 and approaches zero at vanishing temperature 
in consistent with the third thermodynamical law. This rapid change of the entropy is related 
with the Schottky-type specific heat shown in the inset of Fig. [2} 

2. The asymmetric case 

Next we apply the combined method to model A with the asymmetric DW potential 
given by 

U(x) = ^{x^-xlf-d^-xlx^ (28) 
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where d signifies a degree of the asymmetry. Locally-stable minima of the potential locate 
at x — ±x and an unstable maximum is at x u — d (2xq/ Wq) with 

U(±x ) = ± (29) 



U{x u ) 



3 

mw\x\ 



+ 



d 2 xl 



muj. 



o 



AU = U(x ) - U(-x ) 
The asymmetry parameter d is assumed to be given by 

— d r < d < d, 



2d 4 x A 
Adxl 



2V 



(30) 
(31) 

(32) 



for which x u locates at — Xq < x u < Xq. We obtain d c = 1/7 for adopted parameters of 
m = 1.0, u = 1.0 and x = 3.5. In the limit of d — 0, U{x) in Eq. (28) reduces to the 



symmetric DW potential given by Eq. (23) 



Table 1 shows potential values of U(—xq), U (x u ), U(xq) and AU as a function of d. When 
a sign of d is changed, those of U(—Xq), U(xq) and AU are changed, but U (x u ) is unchanged. 
With increasing \d\, \AU\ is gradually increased. Chain and double-chain curves in Fig. [T] 
show U(x) for d = 0.01 and —0.01, respectively, for which a difference of \AU\ is about 30 
% of the potential barrier of \U(x u ) — U(—Xq)\. 



Matrix elements H mn of Eq. ( 17) are not vanishing for pairs of \m— n\ < 4. Eigenvalues for 
n = and 1 are quasi-degenerated for d = 0.0 as mentioned before. This quasi-degeneracy is 
removed with an introduction of d: S (= E\ — Eq) is increased with increasing \d\ as shown in 
Table 1. Eigenvalues e n for d = 0.0 (circles), 0.005 (triangles) and 0.01 (squares) evaluated 
by the energy-matrix diagonalization with N m = 30 are plotted as a function of n in Fig. |5j 
where an increase in 5 with increasing d is clearly realized. 



d 


U(-x Q ) 


U{x u ) 


U(x ) 


AU 


5 


0.0 


0.0 


1.53125 


0.0 


0.0 


0.001943 


±0.001 


^0.02858 


1.53140 


±0.02858 


±0.057166 


0.053123 


±0.002 


^0.05716 


1.53185 


±0.05716 


±0.11433 


0.10619 


±0.003 


^0.08575 


1.53260 


±0.08575 


±0.16184 


0.15927 


±0.004 


=F0. 11433 


1.53337 


±0.11433 


±0.22867 


0.21235 


±0.005 


=F0. 14292 


1.53500 


±0.14292 


±0.28583 


0.26542 


±0.01 


^0.28583 


1.54624 


±0.28583 


±0.57167 


0.53070 
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FIG. 5: (Color online) The n dependence of eigenvalues e n of model A with the asymmetric 
DW potential [Eq. p8] )] for various d: d = 0.0 (circles), 0.005 (triangles) and 0.01 (squares) with 
N m = 30, dashed curves being plotted for guide of the eye. 

Table 1 Potential values at locally-stable minima (±xq), an unstable maximum position 
(x u ), AU [= U(xo)—U(—x )], and the energy gap 5 (= ex— e ) as a function of the asymmetry 
d of model A with the asymmetric potential U(x) [Eq. ( p8| )] (N m = 30). 

Figure [6] shows quantum specific heats calculated with asymmetric potentials for various 
d values [18]. The specific heat for d = 0.0 has the Schottky-like anomaly at very low 
temperature of T ~ 0.001 (see the inset). When a small asymmetry of d = 0.001 (or 0.002) 
is introduced, the position of the Schottky-type peak moves to higher temperature because 
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FIG. 6: (Color online) The temperature dependence of the quantum specific heat C(T) of model 



A with the asymmetric potential [Eq. (28)] for various d: d = 0.0 (bold solid curve), 0.001 (chain 
curve), 0.002 (dashed curve), 0.003 (double-chain curve), 0.004 (dotted curve) and 0.005 (solid 
curve), curves being successively shifted upward by 0.25 for clarity of figure. The inset shows an 
enlarged plot of C(T) with d = 0.0 at < kBT/hujQ < 0.01. Note that the temperature dependence 
of C(T) for a negative d is the same as that for a positive \d\. 

of an increased gap of 5. For d = 0.005, the Schottky-type peak almost disappears and its 
trace is realized as a shoulder at ksT/hujQ ~ 0.1. When we adopt a negative d, AU changes 
its sign but 5 does not (Table 1). The temperature dependence of C(T) for a negative d 
with AU < is the same as that for a positive \d\ with AU > 0. Although the asymmetry 
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has appreciable effects on the specific heat at low temperatures, it has no effects at higher 
temperatures of k B T/hu > 1.0 for adopted asymmetry parameters. 



B. A quadratic DW potential perturbed by Gaussian barrier (model B) 

We will apply our combined method to model B with a symmetric quadratic potential 
perturbed by a Gaussian barrier given by [101 [EE] 



U{x) = U (x) + ae- bx + c, (33) 
Uo{x) = -^—i (34) 



where a and c are parameters. The potential given by Eq. (33) has stable minima at 



x = ±^y\n(2ab/mujQ)/b = ±Xo and an unstable maximum at x = 0.0. For our numerical 
calculations, we assume m = 1.0, Uq = 1.0, a = 9.0, b = 1.0, c = —1.945 and Xq = 1.700, 
which yield U{x ) = 0.0, U(0) = 7.055, U"(x ) = 5.78, and A = 17(0) - U(x ) = 7.055. 
The adopted potential is plotted by solid curve in Fig. [7] where dashed curve expresses the 



harmonic potential given by Eq. (34) 



We have numerically calculated the classical partition function to obtain the classical 
specific heat and entropy. The calculated classical specific heat plotted by chain curve in 
Fig. [8]is not in good agreement with Cho of a harmonic oscillator at kgT ' jhw^ < 10 although 
both reduce to C = ks in the high-temperature limit. The calculated classical entropy will 



be explained shortly (Fig. 10). 



For quantum statistical calculation, we have numerically evaluated eigenvalues by the 



energy-matrix diagonalization. Matrix elements H mn of Eq. (17) are finite for any pair of 
even \m — n\. Then the energy-matrix diagonalization for model B is more time consuming 
than that for model A. Eigenvalues e n calculated for N m = 20 and 30 are plotted in the inset 
of Fig. |9](a). Eigenvalues for n < 20 are almost the same for N m = 20 and 30. Eigenvalues e n 
for n = 0, 1, 2 and 3 are 1.13021, 1.13332, 3.1931 and 3.21918, respectively, with N m = 30. 
eo and e\ are quasi-degenerated with a small gap of 5 = e\ — eo = 0.00311. We note that 
eigenvalues e n of model B in the inset of Fig. |9^a) are similar to but slightly different from 
those of model A in the inset of Fig. [3]^a) . 

Quantum specific heats calculated with the use of eigenvalues e n for < n < N m with 
N m = 20 and 30 are plotted by dashed and chain curves, respectively, in Fig. |9^b). Both 
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FIG. 7: (Color online) The symmetric DW potential of model B [Eq. (33)] (solid curve) and the 
harmonic potential [Eq. (p4|)] (dashed curve). 



results with N m = 20 and 30 are in good agreement each other at ksT/hujo < 2 but 
significantly different at kBT/hu > 4. We assume that extrapolated eigenvalues are given 
by 



I 



n 



for N„ 



1 < n < oo. 



(35) 



The combined eigenvalues are shown by the dashed curve in Fig. |9[a). The quantum specific 
heat calculated with the combined eigenvalues is shown by the solid curve in Fig. [8] [or Fig. 
[9Tb)]. The quantum specific heat has the Schottky-type peak at very low temperature at 
hsT/hido ~ 0.0015 ~ S/2, as shown in the inset of Fig. [8j An increase of the quantum 
specific heat with raising the temperature from zero is slower than that of the harmonic 
oscillator plotted by the dashed curve in Fig. |HJ This is due to the fact that the curvature 
of U(x) at the locally stable point is larger than that of Uq(x): U" {xq) (= 5.78) > t/o'(0) (= 
1.0), which is realized in Fig. [7j A comparison between Fig. [8]and Fig. [2]shows that although 
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FIG. 8: (Color online) Temperature dependences of classical (chain curve) and quantum specific 



heats (solid curve) of model B with the symmetric DW potential [Eq. (33)], dashed curve expressing 
quantum specific heat of a harmonic oscillator (HO). The inset shows an enlarged plot of the 
quantum specific heat at very low temperatures with the Schottky-type anomaly. 

C(T) of model B is similar to that of model A at very low temperatures (ksT/huQ < 0.02), 
they are rather different at higher temperatures (ksT/ Huq > 0.5). In the limit of T — > oo, 
we obtain C(T) = ks in model B while C(T) = (3/4) ks in model A. 

The temperature dependence of the classical and quantum entropies of model B are shown 



by chain and solid curves, respectively, in Fig. 10 , where the quantum entropy of a harmonic 



oscillator is plotted by the dashed curve for a comparison. The inset of Fig. 10 shows the 
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FIG. 9: (Color online) (a) The n dependence of eigenvalues e n of model B obtained by the energy- 
matrix diagonalization for < n < N m (circles) and by an extrapolation given by e' n = (n + 
1/2)Hujq for n > N m (dashed curve), the inset showing eigenvalues obtained by the energy-matrix 
diagonalization with N m = 20 (open circles) and N m = 30 (filled circles), (b) The temperature 
dependence of the quantum specific heat calculated with eigenvalues e n for < n < N m with 
N m = 20 (chain curve) and N m = 30 (dashed curve), the solid curve expressing the result with 
combined eigenvalues shown by the dashed curve in (a). 

quantum entropy at very low temperatures. With raising the temperature from zero, the 
entropy is rapidly developed to 0.69 (~ In 2) at /cgT '/hco® ~ 0.005, which is related with the 
Schottky-type specific heat at very low temperatures shown in the inset of Fig. [8j 
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FIG. 10: (Color online) Temperature dependences of classical (chain curve) and quantum entropies 
(solid curve) of model B, and the quantum entropy of a harmonic oscillator (HO), the inset showing 
an enlarged plot of the quantum entropy at very low temperatures. 

C. An asymmetric DW potential (FUK model) 

The specific heat of a DW system was calculated by FUK [15] with the use of ZOM [16] 
which is explained in the Appendix. FUK adopted an asymmetric DW potential given by 

Ufuk{x) = l:X 2 - Ax 3 + 7X 4 , (36) 

with 

|A| > A m = ^ yfy, (37) 
19 



which has locally stable minima at x = and x s = (3 A/87) [1 + \Jl — I67/9A 2 ], and an 
unstable maximum at x u = (3A/87)[1 — v/I — I67/9A 2 ]. Note that a prefactor of x 2 in 
Ufuk(%) is positive which is required for an application of ZOM [16J, while that of the 



quartic DW potential given by Eq. (|28|) is negative. The DW potential given by Eqs. 

: U FUK (0) = and 



(36) and (37) becomes symmetric with respect to x u with Ufuk(x s 
U F uk{xu) = 1/647 for 



A = A r 



27, X u 



1 

2A? 



X s 



1 
A 



For this symmetric case, a change of a variable x with u = 1/2A leads to 

,2 1)2 



U F uk{x + u) 



57X 



64 7 



which is equivalent to U(x) of model A in Eq. (23). 



A 


UfUK\ x u) 


Ufuk{x s ) = AU 


5 


0.06340 


7.7370 


-0.61721 


0.59848 


0.06335 


7.7613 


-0.41593 


0.41286 


0.06330 


7.7857 


-0.21605 


0.23789 


0.063246 (= A c ) 


7.8125 


0.0 


0.11705 


0.06320 


7.8351 


0.17948 


0.20636 


0.06315 


7.8600 


0.37514 


0.37524 


0.06310 


7.8852 


0.56939 


0.55407 



(38) 



(39) 



Table 2 Potentials values at an unstable position (x u ) and a stable position (x s ), the 
potential difference AU [= Ufuk{x s ) — Ufuk(Q)] an d the energy gap 8 (= e\ — e ) as a 
function of A with 7 = 0.002 for U FU k{x) [Eq. d36b] (N m = 30). 



We have applied our combined method to a DW system with Ufuk(x) given by Eq. (36) 
with necessary modifications. We have chosen potential parameters of 7 = 0.002 and various 



A after FUK (see below). The solid curve in Fig. 11 expresses Ufuk(x) for 7 = 0.002 and 
A = ^/2~7 (= A c ~ 0.063246), which is symmetric with respect to x = 7.9057 (= x u = x s /2) 
with Ufuk(x s ) = Ufuk(Q) = 0.0 and Ufuk(x u ) = 7.8125. When A is varied for a fixed value 
of 7 = 0.002, the potential difference between the two minima, AU = Ufuk(x s ) — Ufuk(0), 
changes: AU > (AU < 0) for A < A c (A > A c ), as shown in Table 2. Ufuk(x) for 
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FIG. 11: (Color online) The DW potential U FUK {x) given by Eq. ([36]) for A = 0.06310 (dashed 
curve), 0.063246 (solid curve) and 0.06340 (chain curve) with 7 = 0.002, the inset showing AU 
[= Ufuk{%s) — Ufuk(0)] as a function of A (xlOO). 



A = 0.06310 and 0.06340 are plotted by dashed and chain curves, respectively, in Fig. [TT 
whose inset expresses AU SIS db function of A. 

By using the energy-matrix diagonalization, we have evaluated eigenvalues of {e„} for 
various A with 7 = 0.002 for < n < N m (= 30). Table 2 shows that the calculated energy 
gap of 5 (= e\ — eo) is minimum for A = A c and it is increased with increasing |A — A c |. 
By using obtained eigenvalues, we have calculated specific heats for various A whose results 



are plotted in Fig. 12 For A = A c ~ 0.063246 (solid curve) with AU = 0, the peak of the 
Schottky-type specific heat locates at T ~ 0.05 ~ 5/2. With decreasing A from A c , both AU 
and S are increased, and then the peak position of the Schottky-type specific heat moves 
upward. For A < 0.06315 the peak of the Schottky specific heat disappears, merging with 
the bump. On the other hand, with increasing A from A c , AU becomes negative as shown in 
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FIG. 12: (Color online) Temperature dependences of the quantum specific heat of a DW system 
with Ufuk{x) [Eq. (36)] for 7 = 0.002 with various A: A = 0.06310 (chain curve), 0.06315 (dashed 
curve), 0.06320 (double chain curve), 0.063246 (solid curve), 0.06330 (bold double-chain curve), 
0.06335 (bold dashed curve) and 0.06340 (bold chain curve), curves being successively shifted 
upward by 0.25 for clarity of figure. 



Table 2. We note in Fig. 12 that temperature dependences of C(T) for A = 0.06330 (bold 
double-chain curve), A = 0.06335 (bold dashed curve) and A = 0.06340 (bold chain curve) 
are nearly the same as those for A = 0.06320 (double chain curve), A = 0.06315 (dashed 
curve) and A = 0.06310 (chain curve), respectively. Thus C(T) of the FUK model is nearly 
symmetric with respect to A = A c . This is similar to the case of model A (Fig. [6]) where 
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C(T) is symmetric with respect to a degree of the asymmetry d. 



IV. DISCUSSION 

A. A comparison between the results of FUK and ours 



The specific heat calculated for the FUK model shown in Fig. 12 has been compared with 
Cfuk(T) reported in Fig. 3 of Ref. [15]. A comparison between the two results shows that 
they are quite different in the following points: (i) The Schottky-type anomaly of Cfuk(T) 
locates at much higher temperature with wider width than ours, (ii) Cfuk{T) has more 
complicated temperature dependence than ours, and (iii) the temperature dependence of 
our specific heat is almost symmetric with respect to A = A c , but Cfuk(T) is not. As for 
the item (i), we suppose that it arises from a neglect of off-diagonal term H od in ZOM [see 



Eq. (A2) in the Appendix]. Among matrix elements of the Hamiltonian, H mn , given by 

H mn = (m\H\n) = (n\H d \n)5 mn + (m\H od \n)(l - S mn ), (40) 



only the first diagonal term is included in ZOM [15] whereas both diagonal and off-diagonal 
terms for < n,m < N m are taken into account in our numerical diagonalization method 



[Eqs. (16) and (17)]. We should note that the off-diagonal term plays an essential role 
in yielding a gap between two stable states in the DW potential. Indeed if off-diagonal 
contributions are neglected in our calculation, we cannot obtain the Schottky-type specific 
heat. 

As for the item (ii), FUK claimed that the calculated, complicated T dependence of 
Cfuk{T) originates from the temperature-dependent energy gap [15]. This arises from the 



fact that optimum parameters of u>, u and v given by Eq. (All ) in ZOM are determined at 



each temperature and then <j>(u,u, v) in Eq. (A10) is temperature dependent in general. It 



is natural that Cfuk{T) is different from the Schottky-type specific heat which is obtained 
for the constant (temperature-independent) energy gap. 

Related to the item (iii), FUK pointed out that Cfuk(T) has a singularity when AU 
changes its sign [15J. Such a result is, however, not realized in our calculation. As mentioned 
before, our A-dependent C(T) is almost symmetric with respect to A = A c . On the contrary, 
the temperature dependence of Cfuk(T) for A > A c is quite different from that for A < A c . 
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FUK [15] reported the specific-heat calculation also for a different set of parameters of 
A = 0.01 and 7 = 0.002 for which Ufuk{x) has a single-minimum structure because they 



do not satisfy the condition given by Eq. (37) (A m = 0.0596 for 7 = 0.002). We have 



realized that Cfuk{T) for these parameters (see Fig. 2 of Ref. |15j ) is in agreement with 
the specific heat obtained by our calculation (related results not shown). It is suggested 
that although ZOM is not applicable to DW systems, it may provide reasonable results for 
single- well systems where off-diagonal contributions are expected to be unimportant. This 
is consistent with the fact that ZOM yields good results for systems with the anharmonic 
potential [TH] and the Morse potential 



B. Triple- well systems 

Finally we will study the triple-well system with the sextic potential 

^) = -^ + X + ^"' (41) 

with parameters a, b and c, whose quasi-exact eigenvalues were investigated in Ref. [5J. 
When we concentrate our attention to low-lying three states with U(—x s ) — U(0) — U(x s ), 
their eigenvalues are approximately given by 

e , e! = e + 5, e 2 = e x + 5, (42) 

where x = and ±x s express locally stable positions and 5 denotes energy gap between 
successive states. At low temperatures, this energy spectrum yields the Schottky specific 
heat whose peak locates at ksT = (3/4) 5. On the other hand, at high temperatures, 
the sextic potential leads to the classical specific heat given by C = (1/2 + l/6)k B = 



(2/3)k B . The sextic potential given by Eq. (41) may be of single-, double- and triple- 
minima type, depending on parameters of a, b and c. It would be worthwhile to investigate 
how thermodynamical properties of the sextic potential system are changed against a change 
of the potential type, whose detailed study is left as our future subject. 

V. CONCLUDING REMARK 

We have calculated specific heats of quantum DW systems with a quartic potential (model 
A), a quadratic potential perturbed by Gaussian barrier (model B) and Ufuk(x) (FUK 
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model) [15J, by using the combined method in which eigenvalues obtained by finite-size 
energy-matrix diagonalization as well as extrapolated ones are included. Specific heat and 
entropy in models A and B with symmetric potentials have the Schottky-type anomaly at 
very low temperatures, which arises from low-lying eigenstates with a small gap due to 
a tunneling through the potential barrier. This is a quantum effect characteristic in DW 
systems, which is sensitive to an asymmetry in DW potentials. In the high-temperature 
limit, specific heats of models A and B reduce to C = (3/4) ks and C = k B , respectively: 
the former is different from Cho — ks of a harmonic oscillator. 

Advantages of our numerical combined method are that the calculations is physically 
transparent and that it yields correct results in both low- and high-temperature regions. We 
have pointed out that the specific heat of DW systems calculated with ZOM [15J is incorrect 
because it neglects off-diagonal contributions which play essential roles for a tunneling in 
DW potentials. Although analytical methods such as the path-integral method (PIM) [T2T 
[T4] and the Gaussian wavepacket method (GWM) [20] have been proposed to obtain the 
partition function of quantum DW systems, they are not suited for calculations of their 
specific heat [21]. The present calculation has clarified a basic problem on the specific heat 
of a DW system expressed by a pedagogical toy model, which is a basis for a study on more 
realistic DW systems, for example, described by system-plus-bath models [22| [23]. 
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Appendix: Zeroth-order operator method 



We will briefly mention ZOM [T5] IT5] in which operators p and q are transformed by 



V 



>\l-\SV - a), q 



(a' + a) + u, [a, a'] 



(Al) 



with parameters u and u, and a denoting creation and annihilation operators, respectively. 



Substituting Eq. (Al) to Hamiltonian given by Eq. (36), we obtain 



H = H d + H od , 
25 



(A2) 



where Ha and H d are diagonal and off-diagonal parts, respectively, with [Hd, n] = and 
h = a^a. Neglecting the off-diagonal term H d, we retain only the diagonal term in ZOM, 

H d \ip n ) ~ E n (u,u)\ip n ), \ip n ) ~ \n), a l a\n)=n\n), (A3) 

where \ip n ) and E n (co,u) denote approximate eigenfunction and eigenvalue, respectively, and 
E n (u,u) are expressed in terms of u, u and n (see Eqs. (6)-(8) in Ref. [15J). The partition 
function expressed by 

oo oo 

m = jV*" - Y,^\e- pH \^n) (A4) 

n=0 n=0 

is transformed to 



Z{$) = (v\e~ R \v), (A5) 
where operators R, k and a state |z/) are given by 

R = /3H + k\nv + \nN(u), N(v) = l-v, (A6) 

oo 

\v) = N{vf/ 2 ^ n ' 2 \^n), (u\v) = l, (A7) 

n=0 

k\i/>k) = *#*), (A8) 

with a parameter v (g [0,1]). By using the Bogoljubov inequality, Z(/3) is approximately 
calculated by 

Z(f3) > e-^SM = e-Mw) = Z (p), (A9) 

where 

oo 

-P<f>(u,u,v) = -pN{y) v n En{u, u) - V - In v - ln(l - u). (A10) 



Variational conditions given by 

d(j)(u>,u,v) d(j)(oj,u,v) d(f)(cj,u,v) 



0, (All) 



<9w du dv 

yield self-consistent equations for optimum values of u, u and u, from which the approximate, 
optimized partition function Z (f3) may be obtained. By using ZOM, FUK calculated the 
specific heat of a DW system with Ufuk(x) [15] . 
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